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ABSTRACT 

Results are presented from a 3-D Pluto general circnlation model (GCM) 
that includes conductive heating and cooling, non-local thermodynamic 
eqnilibrinm (non-LTE) heating by methane at 2.3 and 3.3 microns, non- 
LTE cooling by cooling by methane at 7.6 microns, and LTE CO rotational 
line cooling. The GCM also inclndes a treatment of the snbsnrface temper¬ 
ature and surface-atmosphere mass exchange. An initially 1 m thick layer 
of surface nitrogen frost was assumed such that it was large enough to act as 
a large heat sink (compared with the solar heating term) but small enough that the 
water ice subsurface properties were also significant. Structure was found in all three 
directions of the 3-D wind field (with a maximum magnitude of order 10 m/s in the 
horizontal directions and 10“5 microbar/s in the vertical direction). Prograde jets were 
found at several altitudes. The direction of flow over the poles was found to very with 
altitude. Broad regions of up-welling and down-welling were also found. Predictions 
of vertical temperature profiles are provided for the Alice and REX instruments on 
New Horizons, while predictions of light curves are provided for ground-based stellar 
occultation observations. With this model methane concentrations of 0.2% and 1.0% 
and 8 and 24 microbar surface pressures are distinguishable. For ground-based stel¬ 
lar occultations, a detectable difference exists between light curves with the different 
methane concentrations, but not for different initial global mean surface pressures. 

Key words: methods: numerical - occultations - Kuiper belt objects: individual: 
Pluto - planets and satellites: atmospheres - hydrodynamics 


1 INTRODUCTION 

The modelling of Pluto’s atmosphere began with its discov¬ 
ery in 1988 by the stellar occultation technique (Elliot et al. 
1989). These ranged from mathematically convenient mod¬ 
els, such as that of Elliot & Young (1992), which assumed 
that the temperature dependence with height was a power 
law, to those of Yelle & Lunine (1989), which used conduc¬ 
tion and radiative transfer in the non-local thermodynamic 
equilibrium (non-LTE) regime to describe heating by CH 4 at 
3.3 fim and cooling by CH 4 at 7.6 ftm in order to predict 
temperature with height. This type of radiative-conductive 
model was expanded to an additional CH 4 line at 2.3 fim and 
LTE cooling by CO rotational line emission by Strobel et al. 
(1996), and a troposphere with a user-specified height was 
added by Zalucha et al. (2011b). Zhu et al. (2014) further 
expanded this model to the upper atmosphere of Pluto by 
adding the processes of atmospheric escape, UV heating, and 
adiabatic cooling due to atmospheric expansion. Meanwhile, 
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a technique was developed by Elliot et al. (2003a) to invert 
stellar occultation light curve flux vs. time to temperature 
vs. height. This algorithm required a boundary condition 
at the “top” of the atmosphere (usually ~ 100 km altitude, 
depending on the surface radius assumed), which was then 
used to integrate downward in the atmosphere (closer to the 
occultation midtime) assuming hydrostatic equilibrium, the 
ideal gas law, the linear proportionality of number density 
to refractivity, and the equations of light in the geometric 
limit of optics. 

The general result of 1-D modelling efforts (e.g., Yelle 
& Lunine 1989; Hubbard et al. 1990; Elliot & Young 1992; 
Strobel et al. 1996; Hansen & Paige 1996; Elliot et al. 2003a; 
Zalucha et al. 2011a,b; Young 2012; Young 2013; Zhu et al. 
2014), spectroscopy (e.g., Owen et al. 1993; Trykaet al. 1993; 
Lellouch et al. 2009; Greaves et al. 2011; Lellouch et al. 
2011; Holler et al. 2014), and stellar occultation observa¬ 
tions (e.g., Elliot et al. 1989, 2003b; Sicardy et al. 2003; 
Pasachoff et al. 2005; Elliot et al. 2007; Young et al. 2008; 
Person et al. 2008) was that Pluto’s atmosphere was pri¬ 
marily composed of N 2 , CH 4 , CO, and ethane. Pluto’s sur- 
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face temperature was 40 ± 2 (Tryka et al. 1994). Before the 
flyby of Pluto by NASA’s New Horizon’s spacecraft on 2015 
July 14, evidence was inconclusive about the presence of a 
troposphere and its depth (Stansberry et al. 1994; Zalucha 
et al. 2011b). Preliminary results from New Horizons show 
a shallow troposphere confined to the boundary layer (Stern 
2015). Otherwise, the lower atmosphere (referred to in Pluto 
literature as the stratosphere) is characterized by a sharp 
temperature inversion (temperature increasing with height) 
of ~ 5 K km“^). The temperature becomes mostly isother¬ 
mal at altitudes of ~ 100-600 km. Preliminary results from 
New Horizons have determined the surface radius of Pluto 
to be 1187±4 km and the global mean surface pressure to 
be 10 /ibar (Stern 2015)^. Prior to the flyby, analysis of stel¬ 
lar occultation data and CH 4 spectral lines estimated the 
global mean surface pressure to be between 6 and 24 /rbar 
and changing with time (Lellouch et al. 2009; Zalucha et al. 
2011a; Hansen & Paige 1996; Young 2013). Specifically, the 
pressure underwent an increase, perhaps by as much as a 
factor of two from 1988 to 2003, and remained mostly con¬ 
stant after 2003. Based on modelling studies, it was hypoth¬ 
esized that Pluto’s atmosphere sublimates and condenses in 
vapour pressure equilibrium with surface ice as Pluto orbits 
the Sun. Preliminary results from New Horizons also show 
possible evidence for sublimation (Stern 2015). 

An important aspect of those 1-D models is that they 
ignored the transport of heat by wind and adiabatic heating 
and cooling due to downward and upward motions of air. 
There is no a priori reason to suspect that Pluto’s atmo¬ 
sphere lacks wind (and preliminary results from New Hori¬ 
zons have tentatively identihed aeolian surface features Stern 
2015); thus, modellers have begun investigating wind using 
general circulation models (GCMs) for both Pluto and Tri¬ 
ton (Mueller-Wodarg et al. 2001; Vangvichith & Forget 2011; 
Michaels & Young 2011; Miller et al. 2011; Zalucha & Gulbis 
2012; Zalucha & Michaels 2013; Toigo et al. 2015). GCMs 
solve the Navier-Stokes equations, continuity of mass and 
energy, and an equation of state on a global scale with a 
horizontal resolution of a few degrees and a vertical range of 
several scale heights. They are advantageous because they 
solve the system of equations describing geophysical fluid 
dynamics, but disadvantageous because they are computa¬ 
tionally expensive. 

The list of citations above represent six different efforts 
at transforming terrestrial GCMs to Pluto and/or Triton 
(which have similar atmospheres) using knowledge of other 
solar system bodies with atmospheres where applicable. The 
mystery of Pluto’s wind structure at present cannot be an¬ 
swered observationally, since the effects of wind can only be 
ascertained remotely in certain situations (e.g. cloud move¬ 
ments, ocean roughness, and/or surface dunes). Thus, GCMs 
are extremely important for investigating Pluto’s atmo¬ 
spheric circulation with full global and temporal coverage, 
which in turn affects its 3-D temperature field, surface pres¬ 
sure, and transport of volatiles. Currently there is no com¬ 
mon consensus among these GCMs as to Pluto’s circulation. 
Even with the same input parameters, GCMs may produce 
different results due to differing numerical schemes at the 
fundamental level. In more mature helds such as Earth or 
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Martian GCMs, intercomparison studies are conducted that 
strip GCMs down to their basic components. The field of 
Pluto GCM research is simply too new to have had any of 
these studies performed (though, it appears there is good 
reason to do so in the future). 

The two main instruments on New Horizons that ob¬ 
served Pluto’s lower atmosphere are Alice (Stern et al. 2008) 
and the Radio science Experiment (REX, Tyler et al. 2008). 
Alice is an ultraviolet imaging spectrometer. Alice has two 
modes: one that observed the airglow emission (a collection 
of processes that occur in the upper atmosphere at night) of 
Pluto’s atmosphere and one that performed occultation ex¬ 
periments. The particular configuration of the Pluto GCM 
(PGCM) used in this study cannot predict airglow; how¬ 
ever, the latter capability is of interest here. The first type 
of occultation occurred when Pluto occulted the Sun, such 
that sunlight shined through Pluto’s atmosphere and was 
received by New Horizons. The second type of occultation 
was similar, except Pluto occulted a background star. These 
types of occultations, described in Section 5, allow tempera¬ 
ture to be derived from light intensity based on the attenua¬ 
tion of light due to differential refraction of light through an 
atmosphere. Although the flyby has already occurred, it will 
take some time for the data to be transmitted to Earth and 
subsequently reduced to a form that can be compared with 
physical models. The predictions of this paper are therefore 
relevant. 

REX performed a different type of occultation exper¬ 
iment. REX received radio waves sent from Earth as they 
passed through Pluto’s atmosphere. The observed frequency 
modulation, again due to differential refraction of light 
through an atmosphere, also allowed for the retrieval of tem¬ 
perature with altitude. A key point of REX’s observations 
was that the temperature was measured down to the sur¬ 
face (and the surface radius of Pluto itself was for the hrst 
time well-constrained), something that ground-based stel¬ 
lar occultations cannot do. For Pluto, ground-based stellar 
occultations measure down to about 20 km altitude (e.g. Za¬ 
lucha et al. 2011a). This inability to measure temperature 
to the surface has to do with the fact that the S/N ratio 
near the midtime of the light curve becomes too low, and 
these observations essentially “run out of light.” The fun¬ 
damental difference between UV, solar, or IR occultations 
(as typically performed in ground-based studies) vs. radio 
occultations is the latter provided a temperature measure¬ 
ment near the surface for the first time. 

In this paper, the development of the Zalucha & Gulbis 
(2012) (hereinafter known as ZG12) and Zalucha & Michaels 
(2013) (hereinafter known as ZM13) line of Pluto GCMs 
(PGCMs) is continued. The former was a 2-D PGCM (lat¬ 
itude, height, and time) and the latter was a 3-D PGCM 
(longitude, latitude, height, and time). The main difference 
with adding the longitude direction was that (transverse) 
waves in the longitudinal direction were allowed. This ex¬ 
tra dimension allowed for the formation of thermal tides, 
Rossby waves, Kelvin waves, and gravity (buoyancy) waves. 
These waves are physically real and can transport momen¬ 
tum and heat in all three dimensions as well as dissipate 
artificial noise produced by the model. ZM13 predicted a 
prograde jet at the equator with a maximum magnitude of 
10-12 m s“^. The surface pressure predicted for 1988, when 
compared with observed stellar occultation light curves, was 
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8±4 //bar. Likewise for 2002 and 2006 it was 20±4 //bar, 
and for 2007 it was 16±4 //bar. All years had a greater than 
1% concentration of CH 4 needed to fit the data, which was 
higher than results derived from modelling of spectroscopic 
observations (0.5% Lellouch et al. 2009) and the preliminary 
value from New Horizons of 0.25% (Stern 2015). 

The primary upgrade to the PGCM of ZM13 
was to replace the radiative-conductive model with 
that of ?, which includes a more sophisticated treat¬ 
ment of the CH 4 3.3 and 7.6 //m lines, the addition 
of the CH 4 2.3 //m line, and the effects of CO cool¬ 
ing. Furthermore, a new method of equilibrating the 
atmosphere during the initial spin-up time (devel¬ 
oped independently of Toigo et al. 2015) has been 
implemented. A surface-atmosphere mass exchange 
prescription and subsurface and surface model have 
also been introduced, but have not been extensively 
tested (a subject of future work). In Section 2 the model 
is described, in Section 3 the simulation setup is described, in 
Section 4 simulation results are shown for key variables (tem¬ 
perature and wind) and interpretation of the model results 
is provided, in Section 5 predictions of visible-wavelength 
stellar occultation light curves that ground-based observers 
would have expected to see at the time of the New Horizons 
encounter are made, and in Section 6 predictions of tem¬ 
perature profiles that were be observed by Alice and REX 
(but not yet downloaded to Earth) are presented. The date 
for the ground-based predictions is arbitrary, as at present 
no ground-based stellar occultations were known occur at 
that time. At present, it is not known if Pluto’s atmosphere 
will stay in this state for several decades (given Pluto’s 248 
Earth-year long year and sluggish atmosphere) or if there 
will be a sudden atmospheric collapse (Young 2013; Olkin 
et al. 2015). In Section 8, a discussion is provided about 
speculations future regarding multi-year simulations. 


2 THE MODEL 


The PGCM is based off of the Massachusetts Institute of 
Technology (MIT) atmospheric GCM (Marshall et al. 1997). 
The model solves the primitive equations of geophysical 
fluid dynamics in 3-D using the finite volume method on an 
Arakawa C-grid. These equations are: horizontal momentum 


Du , 

Dt dx 
Dv , 


Fy, 


( 1 ) 

( 2 ) 


where D\ \/Dt = 9[ \/dt -I- V • [ ] is the total derivative, u 
is the velocity (relative to the body’s surface) in the lon¬ 
gitudinal direction, v is the velocity in the latitudinal di¬ 
rection, X is linear distance in the eastward direction, y is 
linear distance in the northward direction, t is time, <I> is 
the geopotential {gz, where g is the gravitational acceler¬ 
ation at the surface and 2 : is height), / = 2Qsm(j> is the 
Coriolis parameter (where Q, is the body’s rotation rate and 
(j) is latitude), and Fx and Fy are external forcings, in this 
case friction and velocity drag; vertical momentum in the 
hydrostatic approximation 


dp 


= -a, 


( 3 ) 


where p is pressure and a is the inverse of density; an equa¬ 
tion of state, in this case the ideal gas law 


pa = RT, 


( 4 ) 


where R is the specific gas constant (the universal gas con¬ 
stant divided by the weighted average of the molecular 
weights of the atmospheric constituents) and T is tempera¬ 
ture; mass equation 


du dv dui 
dx~^ dy~^ dp 


( 5 ) 


where uj is the vertical velocity in pressure coordinates and 
M is the mass source or sink term due to condensation or 
sublimation; and an energy equation 


DT Da 


( 6 ) 


where c„ the specific heat at constant volume and J is an 
external heating and cooling term described below. 

In the MIT GCM, the vertical coordinate used is not ex¬ 
actly pressure as described in equations (l)-(6), but instead 
a modified pressure coordinate referred to as p* (Adcroft 
& Campin 2004), or more commonly known as p in atmo¬ 
spheric literature. The definition of this coordinate is 


P*{x,y,t) =p{x,y,t) 


P°s{x,y) 

Ps{x,y,t). 


( 7 ) 


For a given x, y, and p grid box, the vertical boundaries of 
the box may move up or down, thus shrinking or expanding 
the box. Figure 1 shows the grid at the start of integration. 
Initially, the vertical boundaries of each vertical grid point 
lie on the same pressure level (left panel). The initial surface 
pressure is globally constant. As the model evolves (right 
panel), the boxes vertically shrink or expand depending on 
the atmospheric dynamics. Note that the surface pressure 
Ps changes with time and in the figure is represented by the 
column height. In practice, the fluctuations of the vertical 
boundaries are small (< |0.1|) compared with their origi¬ 
nal position. Table 1 shows the values of the initial pressure 
levels (the case numbers will be defined in Section 3). The 
MIT GCM includes two sets of staggered grids in vertical, 
the pf grid (“f” for faces) and pc grid (“c” for centres). The 
Pf grid has 31 points, with the top point being at zero pres¬ 
sure and the bottom point being the surface pressure. The 
Pc grid points (of which there are 30) are located halfway be¬ 
tween the Pf grid points. The number of vertical grid points 
was chosen such that vertical structure could be sufficiently 
resolved while making the model computationally feasible. 
The reason for the two grids is that for the purposes of tak¬ 
ing derivatives; T, u, v, and 4? are located on the pc grid and 
cj is located on the pf grid. This definition leads us to the 
model’s upper and lower boundary conditions: a; = 0 (air 
cannot flow out of the atmosphere or into the ground). As is 
common in GCMs, more points are located near the surface 
since the dynamics are more intricate (owing to boundary 
layer physics). The model is insensitive to the precise posi¬ 
tion of the points. This practice also makes sense for Pluto 
given the large vertical temperature gradient in the lower 
atmosphere. 

The forcing terms in equations (1) and (2) include two 
effects. The first is surface friction, represented by a simple 
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t=0 t=t^ 




Figure 1. Schematic of p* coordinate system. For the sake of brevity, only the x-axis is shown but the same process applies to the 
horizontal y-axis. At time t = 0 (left panel), pressure values are defined on the vertical grid such that they are the same globally for a 
given level. The initial surface pressure (p^) is also globally constant. As the model evolves to some time ti later (right panel), air has 
moved between the boxes (due to the physical and thermodynamic forces acting on them) and their vertical boundaries have changed. 
The surface pressure ps is the sum of the height of the column, as depicted in the figure. The variations in the vertical boundaries are 
usually small. 


drag law —k^u and —kvV, where 

kjj = kf max (0, A — ^ . (8) 

\ Ps ~ Pb J 

Here pb is the pressure top of the frictional layer and Ps is the 
surface pressure. The variable fc/ is an an empirical constant 
and acts as a type of bulk coefficient of friction. In more in¬ 
volved treatments of the boundary layer, this forcing term is 
replaced with yet more complex empirical treatments of the 
surface-air interaction, because the scale of surface features 
is usually much smaller than the GCM horizontal grid. The 
only bodies for which we have empirical data to estimate this 
constant are Earth and Mars. For Earth, this parametriza- 
tion (equation 8) was first used in the simple Earth model of 
Held & Suarez (1994) with a value of an inverse Earth day. 
Similar values have been used for Mars (e.g. Lewis et al. 
1996; Nayvelt et al. 1997)). Since there there is no empirical 
data for Pluto, nor to the author’s knowledge has there been 
a published study on the effect of changing surface pressure 
on fcy, it is left as one inverse Earth day here. Furthermore 
Pb/Ps = 0.7 is allowed in accordance with the original Earth 
(and Mars) model, as we have no data to constrain it for 
Pluto. Ideally, sensitivity tests should be performed on the 
effect of the value of fc/ and other similarly unknown pa¬ 
rameters; however, due to the large computational expense 
of Pluto simulations, it will be left to future work. Sensitivity 
tests are performed on other major quantities in Section 4.2. 

The second forcing effect in the horizontal momentum 
forcing term is a ’’frictional” term at the model top. In a 
model atmosphere, a rigid lid (i.e., a;=0) is applied at the 
model top that does not exist in the real world. This model 
top can cause spurious wave reflections off the lid that are 


non-physical. In the real atmosphere, upward propagating 
gravity waves break in the upper atmosphere, depositing 
momentum and creating a drag on the flow. The proper¬ 
ties of these waves can vary with time and space. Because 
we have practically no data on this process for Pluto (but 
see Person et al. 2008; Hubbard et al. 2009), it is difficult to 
predict theoretically, and waves are smaller in scale than the 
grid boxes, model simulations were performed with a con¬ 
stant drag coefficient set at different orders of magnitude and 
spanning a different number of top levels, a common prac¬ 
tice in simple GCMs. It was tested as to which configuration 
does not damp the upper atmosphere so much that it acts 
like part of the rigid lid, while preventing wave reflections. 
In general, the coefficient’s magnitude must increase with 
height. It was found that the best value for the coefficients 
for drag in the top three model levels that are from the top 
downwards: 0.0001013799, 0.0000337933, 0.0000112644 s“^ 
In addition, an 8th-order Shapiro filter is included in the 
momentum forcing term to damp subgridscale, numerical 
noise. 

Surface topography is assumed to be flat. It has been 
found from New Horizons data that mountains locally rise 2 
to 3 km above the surrounding orography. It is also hypoth¬ 
esized that the informally named Sputnik Planum is a basin 
filled with CO ice bordered by mostly higher terrain. As the 
height of the mountains is much less than the scale height of 
50 km, they can be assumed to fall under the domain of the 
surface friction term in the PGCM. It would be interesting 
to see how these mountains affect the flow at low-levels in a 
regional (i.e. mesoscale) simulation. 

The energy forcing term in Eq. 6, i.e., atmospheric 
heating and cooling terms, were taken from Strobel et al. 
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Table 1. Pressure values for vertical grids used in the PGCM 


Level 

Cases 1 and 3 
(/i.bar) 

Cases 2 and 4 
(/^bar) 

1 

8.000000 

24.000000 

2 

7.860086 

23.580258 

3 

7.651933 

22.955798 

4 

7.351503 

22.054509 

5 

6.924851 

21.414532 

6 

6.333587 

19.887658 

7 

5.543393 

17.815471 

8 

4.583648 

15.190562 

9 

3.693600 

12.415871 

10 

2.927441 

9.0931562 

11 

2.291569 

7.828516 

12 

1.778336 

6.105858 

13 

1.371895 

4.715346 

14 

1.054218 

3.639138 

15 

0.807902 

2.793179 

16 

0.617827 

2.138594 

17 

0.471519 

1.634020 

18 

0.358990 

1.245765 

19 

0.272382 

0.947059 

20 

0.205579 

0.716942 

21 

0.153848 

0.539140 

22 

0.113559 

0.401111 

23 

0.081964 

0.293284 

24 

0.057051 

0.208522 

25 

0.037467 

0.141777 

26 

0.024502 

0.089876 

27 

0.016861 

0.051204 

28 

0.004944 

0.024945 

29 

0.001969 

0.009769 

30 

0.000317 

0.002829 

31 

0.000000 

0.000000 


(1996). The 1-D heat equation is 


Cpp{r, t) 


dTGt) 

dt 


J__9 f 2j^ dT{r,t) \ 
r 2 dr\ dr ) 


tinetG^ t'j ^ 


( 9 ) 


where Cp is the specific heat at constant pressure, p{r, t) is 
the air density, K is the thermal conductivity, and R„et{r,t) 
is the heating rate. Note that the terms in Eq. 9 are the 
external heating terms, J, in Eq. 6 . For Pluto a primarily 
N 2 atmosphere with smaller amounts of CH 4 and CO was 
assumed. For an N 2 , CH 4 , and CO mixture. 


K{r,t) = KoT{r,t)°‘ (10) 

where Ko = 5.63 x 10“® J m“^ s“^ q, = 

1.12 (Hubbard et al. 1990). The most mathematically com¬ 
plex part of the model is the specification of Rnet{r, t), which 
is the sum of the non-LTE heating rate by solar near-IR 
absorption in the CH 4 2.3 and 3.3 pm vibrational bands, 
the non-LTE cooling rate due to the CH 4 7.6 pm vibra¬ 
tional band, and the LTE cooling rate by CO rotational line 
emission^. In the ZG12 specification of the model, the disk- 


^ New Horizons found evidence of C 2 H 2 , which is currently being 
added to the Rnet{r, t) term as a non-LTE cooling line at 13.7 pm. 
Moreover, the effects of CH 4 non-LTE heating at 1.67 pm have 
also been identified as a key heating source and is also being 
added to the 1-D heating and cooling prescription (X. Zhu, private 
communication). 


averaged insolation was been replaced with a fully longitu¬ 
dinally, latitudinally, and seasonally varying insolation. Fi¬ 
nally, a stratosphere only configuration (i.e., no troposphere 
is present as in Zalucha et al. 2011b) was assumed, which is 
broadly consistent with the findings of New Horizons (Stern 
2015). 

The PGCM is capable of undergoing both mass and en¬ 
ergy changes due to freezing and sublimation by the primary 
atmospheric constituent (as (J-phase N 2 ice), represented by 
M in equation (5). If the atmospheric or surface temperature 
falls below the freezing point, the amount of heat needed to 
return it to the freezing point is calculated, and this heat 
is provided by latent heating due to mass freezing out of 
the atmosphere. This mass is instantaneously removed from 
the atmosphere and added to the surface ice inventory. In 
the case of a surface temperature below freezing, the mass 
is taken out of the lowest model layer and added to the 
surface ice inventory. Here the heat exchange is implicit, 
and treated separately from the radiative-conductive heat¬ 
ing equation (9), following the relation 

Lm = pAAzCp{Tfs - Ts) (11) 

where L — 253000 J kg~^ is the latent heat of sublimation, 
m is the amount of mass accumulated, A is the horizontal 
area of the grid box, Az is the height of the lowest atmo¬ 
spheric layer, and Ts is the surface temperature (here below 
the frost point at the surface, T/s). What the above equa¬ 
tion states is that enough latent heat will be supplied to 
bring the surface temperature back up to the frost point. If 
the surface is ice-covered and the temperature is predicted 
to be greater than the freezing point, the excess energy is 
used to sublimate an amount of mass such that the surface 
temperature returns back down to the freezing point. The 
sublimated mass is added to the lowest model layer. 

There is an additional subtlety to the sublimation and 
deposition at the surface. As mass is subtracted from the at¬ 
mosphere, the frost temperature decreases, and the amount 
of frost needed to accumulate on the surface to bring the 
surface temperature back to the frost point decreases. Thus 
a feedback develops. For a full treatment of this effect, an 
iterative procedure is needed. However, it has been deter¬ 
mined that the error in ignoring this process is insignificant. 

The PGCM of ZM13 assumed that the surface tem¬ 
perature was fixed at the N 2 surface temperature, i.e. that 
there was an infinite reservoir that could absorb heat at the 
surface. Here the surface and subsurface temperatures are 
allowed to vary. The findings of Stern (2015) have recently 
shown that Pluto’s surface ice distribution (N 2 , CH 4 , and 
CO) is very non-uniform. At the time of the simulations in 
this work, it was not known that this was the case, and an 
entirely N 2 ice surface was assumed based on ground-based 
observations (Owen et al. 1993). It has been assumed that 
the thickness of the surface N 2 ice is a “surface veneer”, which 
was also concluded by Stern (2015). A value of 1 m was cho¬ 
sen for the surface ice thickness. Appendix A describes the 
equations of the surface and subsurface model in detail. 

In the horizontal, a cubed-sphere grid (Adcroft et al. 
2004) with 32 x 32 points per cube face is used, equiva¬ 
lent to a grid spacing of 2.8° at the equator. Compared to 
the more common cylindrical projection grid (i.e., a latitude 
and longitude grid), this type of horizontal grid eliminates 
singularities at the poles that force northward and south- 
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ward winds to zero and removes the requirement for artificial 
Fourier filtering in the high latitudes (in order to maintain 
a practical time step). Lastly, a surface radius of 1180 km 
is assumed (Zalucha et al. 2011a). The effect on the model 
results from using this value compared with the the recently 
measured New Horizons value of 1187±4 km is quite in- 
signihcant. The surface radius comes into play in insolation 
terms, gravitational acceleration terms, and the calculation 
of stellar occultation light curves. A difference by as much 
50 km (as would be the case if the surface radius were the 
value of the lowest value in the literature, 1137±7 Tholen & 
Buie 1990; Reinsch et al. 1994) is negligible compared with 
Pluto’s distance from the Sun (tens of AU). Likewise, com¬ 
pared with Pluto’s radius, an error of 50 km is small and 
would also not significantly impact the gravitational accel¬ 
eration at the surface. However, as shown by (Zalucha et al. 
2011a), an error of 50 km would be significant when calcu¬ 
lating stellar occultation light curves. The present difference 
of ~7 km has a slight effect on the calculation of stellar 
occultation light curves. 


3 MODEL SETUP 


For this work, the following Pluto pole convention has been 
adopted: the north pole is in the same hemisphere as the 
Sun’s north pole and longitude increases eastward (see Zan- 
gari 2015, for a detailed discussion of Pluto pole conven¬ 
tions). In that sense, Pluto rotates in retrograde compared 
to Earth and most Solar System planets. As such, prograde 
flow will be westward and plotted as positive values. Other 
directions remain unchanged. Note this convention is the op¬ 
posite of both the New Horizons team convention and the 
current lAU convention. From the years 1985 to 2015, the 
subsolar latitude (or in other words Pluto’s ecliptic longi¬ 
tude) changes in accordance with Pluto’s orbit for a given 
date. Thus, seasonal changes in insolation are accounted for 
in the model. On this date at the time of the New Hori¬ 
zons flyby, the subsolar longitude was 32° and the subsolar 
latitude was 38° . 

The atmosphere of Pluto is thermally “sluggish” in that 
the radiative-conductive time constant is long, or in other 
words responds slowly to changes in radiative-conductive 
forcing (Strobel et al. 1996; Young et al. 2008). The time- 
scale for this forcing r^c may be calculated as (e.g., Showman 
et al. 2008) 


'7”rc 


5T 


CpP 

dF/dz ’ 


( 12 ) 


where F is the net vertical radiative-conductive flux and ST 
represents a Gaussian perturbation to a radiative-conductive 
temperature profile at some altitude. While the choice of 
ST is somewhat arbitrary, dF/dz scales accordingly, and so 
by examining a wide range of perturbation amplitudes and 
thicknesses, it is usually possible to gain at least an order of 
magnitude estimate, if not better (fractions of Earth days). 
The time-scale for Pluto in the lower atmosphere is of order 
of a few weeks, while the time-scale in the upper atmosphere 
can be a significant fraction of an Earth year. 

For Mars and Earth GCMs, where the radiative time 
constant is two to four Earth days and tens of days, re¬ 
spectively, conduction is ignored, and the initial tempera¬ 
ture structure may be set to any value as long as it is not 


too pathological, and the model will respond and equilibrate 
over a reasonable amount of wall clock time. Eor Pluto’s at¬ 
mosphere, simulations starting with global isothermal tem¬ 
peratures of 80 K and globally quiescent winds show at¬ 
mospheric temperature changes of no more than 1 K over 
40 Earth years of integration time. Since a nearly constant 
temperature atmosphere is far from the radiative-conductive 
equilibrium solution and observations, it suggests that at 
least that these initial conditions are far from the spun-up 
state. 

To initialize this three-dimensional model from this 
isothermal state would be prohibitively expensive from a 
computational standpoint. The radiative-conductive equi¬ 
librium (i.e. steady state, quiescent wind) is then a logical 
choice for an initial condition as it may be quickly calcu¬ 
lated and is assumed to be near the fully spun-up, quasi¬ 
equilibrium state. However the temperature gradients in this 
configuration are so strong that a prohibitively small time 
step is required in the initial stages for stability. (GCMs in 
typical usage do not employ variable time stepping meth¬ 
ods, since this practice can result in different amounts of 
solar energy being received from one day to the next.) A 
hybrid method to spin-up the model has been developed as 
follows. Note that this method is conceptually the same as 
that of Toigo et al. (2015); however, their paper did not ap¬ 
pear in press until long after the simulations in this paper 
were complete, so the method was independently arrived at 
by both groups. 

First, the initial temperature was set to a constant 
value globally with quiescent winds. The surface and sub¬ 
surface temperatures were held constant and the surface- 
atmosphere mass exchange scheme was disabled. Then, 
the technique of Newtonian relaxation was used to thermally 
force the model in place of the direct heating and cooling 
rates from the Zhu et al. (2014) model. Newtonian relax¬ 
ation is a method by which the forcing term in the energy 
equation goes as 

dT (T-Te,) 

dt^ r ’ 

where Teg is the radiative-conductive equilibrium temper¬ 
ature and r is the radiative-conductive time constant. A 
value of T = 30 Earth days has been found to require a 
weak enough forcing while affording a minimal amount of 
wall clock time. The model was run in this configuration for 
90 model Earth days. Then, using the MIT GCM’s restart 
feature, the model was restarted, now with the heating and 
cooling rates specihed directly from Zhu et al. (2014) (i.e., 
not Newtonian relaxation), with the surface and subsurface 
evolution, and the surface-atmosphere mass exchange 
scheme. From then the model was run from the year 1985 to 
2015, i.e., spanning the observational epoch of stellar occul¬ 
tation observations and NASA’s New Horizons Pluto flyby. 

As also recognized by Toigo et al. (2015), current com¬ 
puting power is not enough to equilibrate a PGCM with 
both an atmospheric and surface component. 2-D surface 
models, such as that of Hansen & Paige (1996) have shown 
that just the subsurface/surface alone (with a 1-layer at¬ 
mosphere) takes up to four Pluto years (about one thou¬ 
sand Earth years) of simulation time to equilibrate. ZM13 
showed that vertical velocity in the atmosphere was very 
weak, most likely due to the large increase in temperature 
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Table 2. Surface parameters of PGCM simulations 


Case 

number 

CH 4 

concentration (%) 

Surface 

pressure (/ibar) 

1 

1 

8 

2 

1 

24 

3 

0.2 

8 

4 

0.2 

24 


with height in the stratosphere, which makes the atmosphere 
extremely stable against vertical motions. It was thus rea¬ 
soned that any flow that arose from pressure differences due 
to sublimating or condensing surface ice would be confined 
to the lowest atmospheric layer and would thus be decou¬ 
pled from the atmosphere aloft. Without PGCM simulations 
results that started with an arbitrary N2 frost distribution 
and spanned many Pluto years, it cannot be known whether 
or not the surface ice distribution is accurate, and thus it 
will not be focused on for the remainder of this paper. Toigo 
et al. (2015) on the other hand used a different approach; 
they ran a 2-D surface model, similar to as that of Hansen & 
Paige (1996), for many Pluto years, and used that as their 
initial surface condition. Again, without PGCM simulations 
initialized alongside a surface model spanning many years, 
it cannot be known if this method gives accurate results ei¬ 
ther. Similarly, the computational expense limits extensive 
sensitivity tests of the not well-constrained parameters, as 
is the case in many planetary GCMs. Such tests will have to 
wait until a future paper. 

In this study, four different plausible situations are 
considered based on analysis of stellar occultation light 
curves (e.g., Zalucha et al. 2011a), long-term surface mod¬ 
els (e.g., Hansen & Paige 1996; Young 2013), spectroscopic 
data (e.g., Lellouch et al. 2009), and direct imaging (e.g., 
Buie et al. 2010). Table 2 shows the parameters of each ex¬ 
periment (each simulation will be referred to by its simula¬ 
tion number hereinafter). At the time the simulations were 
carried out, the values of global mean surface pressure and 
CH 4 concentration spanned values predicted from ground- 
based results. Case 3 will be the primary topic of discussion 
for the remainder of this paper, because it is the closest the 
values observed by New Horizons (Stern 2015). The results 
of the other cases are still relevant, as they may represent the 
state of Pluto’s atmosphere in the recent past (or perhaps 
near future), but can be found in Appendix B. 

4 PGCM RESULTS AND DISCUSSION 
4.1 Temperature 

The atmospheric temperature shows little longitudinal (i.e. 
diurnal) variation. This property is to be expected since 
Young et al. (2008) and ZM13 found that the time-scale 
for radiative-conductive forcing is much longer than a Pluto 
day. In Fig. 2 the longitudinally averaged temperature as a 
function of altitude and latitude are shown for 2015 July 14. 
For each longitudinally averaged panel, the altitude values 
are calculated by assuming a globally averaged temperature 
and pressure at each level. As also found by ZG12 and ZM13, 
the latitudinal temperature gradient is small. In Cases 1 and 


2, the temperature is warmer than Cases 3 and 4 because 
the bulk effect of a higher CH4 concentration. 

The effect of initial surface pressure on atmospheric 
temperature is somewhat less straightforward. Comparing 
Cases 1 and 2, it would seem that lower initial surface pres¬ 
sure corresponds to higher temperature at a given level. Be¬ 
cause CH4 is represented in the atmosphere by a percentage 
of the total atmospheric mass, a higher initial surface pres¬ 
sure, i.e. more massive atmosphere, means a larger amount 
of CH4 and in turn a warmer atmosphere. The CH4 concen¬ 
tration is the same but the partial pressure (i.e. the num¬ 
ber of molecules) of CH4 has increased, thus increasing the 
heating rate. Yet, the opposite appears to be true. This con¬ 
tradiction is an artefact of the way temperature has been 
plotted. Had Cases 1 and 2 been plotted with pressure as 
the ordinate, it would be apparent that higher initial sur¬ 
face pressure corresponds to warmer temperatures at all al¬ 
titudes. At altitudes of ~250 to 500 km, the temperature is 
almost isothermal, except for subtle effect of heating by the 
CH4 2.3 and 3.3 /rm lines and cooling by the 7.6 /rm line. 
The temperature profile “wobbles” by about 10 km in this 
region, and the altitudes of local maxima and minima occur 
at different locations (see Fig. 8 of Zalucha et al. 2011a, for 
an example) depending on the initial surface pressure. 

4.2 Longitudinal Wind 

In Fig. 3 the longitudinally averaged longitudinal wind 
(westward positive and prograde) is shown as a function 
of altitude and latitude for 2015 July 14. In all cases, the 
wind in the lowest 10 km is light (less than 2 m s“^). This 
bulk measurement does not rule out small-scale phenomena 
such as plumes, dust devils, or flows from local topographic 
forcing; it is simply that the PGCM has a resolution of 50- 
60 km, which is too coarse to resolve such flows. Just above 
10 km altitude the winds are westward at all latitudes un¬ 
til about 60 km altitude where they change sign near the 
equator (except in the bottom left panel of Fig. 3, where 
the sign change occurs at 100 km altitude). This equatorial 
region of eastward winds extends from about ±30°latitude. 
Note that this is approximately the width of Earth’s tropics, 
but we believe this delineation is purely coincidental, since 
Pluto has a much different rotation rate, size, and horizontal 
temperature gradient compared with Earth. The eastward 
winds are less than the westward winds, and the former do 
not exceed 4 m s”^ in magnitude at any altitude. 

In all cases, the winds are westward in the southern (i.e., 
winter) hemisphere south of —30°latitude, with the excep¬ 
tion of a narrow range of eastward winds around 150 km al¬ 
titude for 8 /ibar initial surface pressure and 0.2% CH4 con¬ 
centration (Case 3). The maximum westward wind occurs 
between —75 and —70° latitude at an altitude of 50-70 km 
for CH4 concentrations of 0.2%. The maximum wind speed 
is of order 10 m s“^. The winds in the northern (i.e. sum¬ 
mer) hemisphere north of 30° latitude are mostly westward 
as well, but slightly weaker. In all cases there is a region 
of weak (~2 m s“^) eastward wind aloft between —20 and 
20°latitude that varies in its altitudinal extent. The alti¬ 
tudes and latitudes of local maxima and minima are mostly 
symmetric with the southern hemisphere. 

Note that the structure of longitudinal winds does not 
at all resemble that of ZM13. While the maximum westward 
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Figure 2. Longitudinally averaged temperature (K) for 2015 July 14. Top left: Case 1, top right: Case 2, bottom left: Case 3, bottom 
right: Case 4. The horizontal temperature gradient is weak compared with other Solar System bodies with atmospheres. Cases 1 and 2 
are the warmest temperature due to the higher CH 4 concentration (which has a net heating effect on the atmosphere) and the fact that 
altitude is plotted on the ordinate rather than pressure. Likewise, for Cases 3 and 4, the temperature is colder for a given altitude than 
in Cases 1 and 2 because of the lower CH 4 concentration. 


and eastward winds (12 m s“^and 2 m s“^, respectively) are 
of the same order of magnitude, ZM13 found a westward 
maximum at the equator and 175 km altitude. Between 50 
and 175 km altitude there was a region of weak eastward 
winds northward of —30° latitude, and weak westward winds 
south of this latitude. Below 50 km altitudes the winds were 
light and lacked structure with regards to direction. The only 
similarity is that the maximum of the westward winds is of 
the same order of magnitude of the present study. The differ¬ 
ences between the model configurations are that ZM13 used 
the radiative-conductive scheme of Yelle & Lunine (1989) 
(no CO or 2.3 fim CH 4 line), the drag at the model top that 
was several orders of magnitude stronger and no surface- 
atmosphere mass exchange scheme was used. The ef¬ 
fects of these processes on the circulation are examined as 
follows. 

Figure 4 shows the longitudinally averaged westward 
wind for Case 3, each with one of the model elements dis¬ 
abled, as well as an additional case run at Northern Hemi¬ 
sphere winter solstice and the full model results (hereinafter 
known as operational results) shown above for easy compar¬ 
ison. The scenarios considered are: (a) no atmospheric CO 
cooling in the radiative-conductive scheme, (b) no contribu¬ 


tion from the 2.3 fim CH 4 line, (c) an upper atmospheric 
drag coefficient that is an order of magnitude higher than 
the operational runs and encompasses the top five layers in¬ 
stead of the top 3, (d) no mass exchange between sur¬ 
face and atmosphere, and (e) initialization at Northern 
Hemisphere winter solstice (integrated for the same amount 
of time as the operational runs). 

The most striking difference is scenario (e), in which the 
winds are the opposite direction (i.e. retrograde or eastward) 
although of the same order of magnitude as the operational 
runs, maximum at the equator, and approximately symmet¬ 
ric at the equator. Note the low-level area of westward winds 
close to the surface remains. Until a follow up study is done 
with full Pluto-year simulations, it is not possible to know if 
this retrograde wind pattern is characteristic for this time or 
if there is a dependency on the model start time. The latter 
is a property that will certainly be investigated in the 
future. 

Another difference from the operational run is appar¬ 
ent in scenario (c). As expected, instead of the winds in 
the upper part of the domain going static around 300 km 
altitude, they now are damped down to altitudes between 
200-250 km. The circulation maxima near the poles are com- 
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Figure 3. Longitudinally averaged westward (prograde) wind (m s“^) for 2015 July 14. Top left: Case 1, top right: Case 2, bottom 
left: Case 3, bottom right: Case 4. Positive values are westward (prograde). For Case 1, the wind is strongest in the westward flow 
at —75° latitude between 50 and 100 km altitude. Around 50 km altitude, the winds are westward at all latitudes, while above this 
layer there is an area of weak eastward winds at the equator. For Case 2, the structure of the winds is the same as Case 1, except the 
maximum westward winds are weaker and the maximum eastward winds are stronger. Case differs 3 the most from the other cases in 
that between roughly 200 and 250 km altitude, the winds are eastward everywhere. The maximum eastward wind is also located between 
these altitudes, higher than the other cases. Furthermore, the westward wind maximum is the strongest of the four cases. For Case 4, 
the wind structure is similar to Case 1, except the westward maximum is slightly stronger, and the eastward maximum is located at a 
higher altitude. 


pacted; however note that the lowest level flow just above 
the surface is not impacted. Unexpectedly, no drastic change 
in the circulation occurred due to purported waves reflect¬ 
ing off of the rigid lid caused by the substantially increased 
damping, except for in Case 3 where a thin band of eastward 
wind appears at about 225 km altitude over the equator and 
midlatitudes. 

Scenario (b) is also clearly different in that the circula¬ 
tion splits into two distinct polar, westward jets with little 
to no flow at the equator and the disappearance of the near¬ 
surface westward flow that is present at all latitudes in the 
operational run. Briefly note that in Scenario (a), where the 
CO cooling has been disabled, the circulation hardly dif¬ 
fers from the operational run (except for the weakening of 
westward winds at the equator at 300 km altitudes). Thus, 
Scenario (b) essentially recovers the case of Yelle & Lu- 
nine (1989). These authors had found that for the radiative- 
conductive equilibrium temperature, the 2.3 fim CH 4 line 


was unimportant. Later, Strobel et al. (1996) argued that 
is was important to the radiative-conductive temperature. 
A further difference between these papers is the former as¬ 
sumes optically thin radiative transfer for near and thermal 
IR in the gray approximation for CH 4 bands and neglect of 
vibrational energy transfer from stretch modes to bending 
modes. In any case, the results of Scenario (b) show that 
assumptions in the radiative-conductive scheme propagate 
into significant differences in the atmospheric circulation. 

Scenario (b) is quite curious because except for Case 
3, the polar structure is mainly preserved, but the westward 
winds intensify over the equator at altitudes between 200 
and 250 km. Again Case 3 stands out because an eastward 
area of wind develops over the equator. Why the prefer¬ 
ence for Case 3 to develop these eastward jets is not easily 
inferred. The near-surface westward winds remain, which 
calls into question whether these are related to the conden- 
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Figure 4. Longitudinally averaged westward (prograde) wind (m s“^) for Case 3 on 2015 July 14 Positive values are westward (prograde). 
The panels are (a) no atmospheric CO cooling in the radiative-conductive scheme, (b) no contribution from the 2.3 /rm CH 4 line, (c) an 
upper atmospheric drag coefficient that is an order of magnitude higher than the operational runs and encompasses the top five layers 
instead of the top 3, (d) no mass exchange between surface and atmosphere, and (e) initialization at Northern Hemisphere winter 
solstice. The panel with the adjacent colorbar is the same as Fig. 3, bottom left panel for comparison. 


sation flow. Figures B1-B3 show the corresponding results 
for Cases 1, 2, and 4. 

Figure 5 show horizontals cross-sections of longitndi- 
nal wind at key altitudes for Case 3 on 2015 July 14. Near 
the surface at 2.2 km altitude (upper left panels), the wind 
is weak, but there is a distinct wave number -1 pattern be¬ 
tween westward and eastward directions. At the next alti¬ 
tude shown (32 km) there is a strong westward jet that is 
of constant strength with longitude. Elsewhere (except very 
near the north pole), the winds are also westward, but there 
is a local minima near 270° longitude and a local maxima 
near 90°. At 98 km altitude, the winds are weakly east¬ 
ward between —30° and 30° latitude (except for an area of 
westward winds in the 24 jrbar, 1% CH 4 case) with a slight 
variation with longitude. The winds are relatively strongly 
westward near the poles. Finally, at 217 km altitude, there is 
a wave number -1 structure between —60° and 60° latitude, 
with a local maximum of eastward winds at 270° longitude 


and a local maximum of westward winds at 90° longitude. 
Poleward of these latitudes, the location of the eastward and 
westward maxima are shifted by 180°. Figures B4-B6 show 
the corresponding results for Cases 1,2, and 4. 


4.3 Latitudinal and Vertical Wind 

In the longitudinal average, the latitudinal and vertical wind 
do not show any structure. However, when plotted in the 
xy-plane, at the same altitudes as Fig. 5, each show distinct 
variations. Figure 6 shows the latitudinal wind for Case 3 
(Figs. B7, B 8 , and B9 show the corresponding results for 
Cases 1, 2, and 4). At the equator, the latitudinal flow is 
nearly motionless. At low altitudes, again the flow is weak for 
the same reason it was in the longitudinal direction. There 
are definitive areas of northward (positive) and southward 
(negative) flow near the poles (poleward of 60° latitude). 
For the uppermost three altitudes shown, between approxi- 
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Longitude (degrees) Longitude (degrees) 

Figure 5. Instantaneous westward (prograde) wind (m s“^) for Case 3 on 2015 July 14. Positive values are westward (prograde). On 
this date at the time of the New Horizons flyby, the subsolar longitude was 32° and the subsolar latitude was 38° (in the pole convention 
used in this work, see text). In the top left panel (2.2 km altitude), the magnitude of the maximum wind is of order 1 m s“^, and the 
eastern (western) hemisphere is eastward (westward). In the top right panel (32 km altitude), there is a strong flow (jet) of westward 
wind near the south pole, while elsewhere there is a local maxima of westward wind in the eastern hemisphere and a local minima in the 
western hemisphere. In the lower left panel (98 km altitude) and lower right panel (248 km altitude), again there is a jet of westward 
wind near the south pole. There is a weaker flow of westward winds near the north pole (and a localized area of eastward wind very close 
to the north and south poles). At 217 km altitude there is an area of weak westward flow in the eastern hemisphere. 


mately 90° and 270° longitude the flow is southward at the 
north pole and northward at the south pole; i.e. away from 
the pole. At other longitudes, the flow is towards the pole. 
For a given longitude, the flow direction is opposite that 
at the opposite pole (e.g., where there is northward flow at 
the south pole, there is southward flow at the north pole 
and vice-versa). Remembering that we are on a 3-D sphere, 
these figures suggest that air is flowing over the poles. 

When all 30 levels of the model are examined (not 
shown), the areas of northward and southward flow at the 
poles shift eastward with altitude. Thus, the air always flows 
over the poles (except near the surface), but the meridian it 
follows shifts eastward with height. In Fig. 5, the altitudes 
chosen all show northward and southward flow confined to 
the same longitudes. Looking at Case 1 (Fig. B4, the top 
two altitudes plotted (bottom two panels) show the areas of 
maximum northward and southward flow 180° out of phase 
in longitude. That is, in the lower left panel between lon¬ 
gitudes 90° and 225° the flow is southward at the south 
pole, but in the lower right panel the flow is northward at 


these same longitudes. The reasons that the areas of 
maximum north and south flow shift eastward with 
height is not yet understood. 


Figure 7 shows the vertical wind in pressure coor¬ 
dinates for Case 3. In this coordinate system, positive 
is downwards, negative is upwards, and the units are in 
10“®/ibar s“^. There is a clear area of up-welling between 
90° and 270° longitude, with down-welling at other longi¬ 
tudes, except at the highest altitude shown, where there is 
little vertical motion. There is also fine-scale structure that 
is almost as large in magnitude as the maxima and minima. 
What we may conclude from the meridional and vertical 
wind directions is that the flow on Pluto is not like other 
bodies in the solar system that show overturning circula¬ 
tions in the longitudinal mean, i.e. Hadley cells. Future work 
explaining the unusual atmospheric dynamics is needed. Fig¬ 
ures B10-B12 show the corresponding results for Cases 1, 2, 
and 4. 
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Figure 6. Instantaneous northward wind (m s“'^) for Case 3 on 2015 July 14. Positive values are northward. On this date at the time 
of the New Horizons flyby, the subsolar longitude was 32° and the subsolar latitude was 38° (in the pole convention used in this work, 
see text). Top left panel: 2.2 km altitude, top right panel: 32 km altitude, bottom left panel: 98 km altitude, bottom right panel: 217 km 
altitude. On this date, the subsolar point is located at 32° longitude. Except at the lowest level, air appears to be flowing over the poles 
along the terminators. The magnitude of the maximum flow is about 10m s"'^. 


5 MODEL STELLAR OCCULTATION LIGHT 

CURVES 

Stellar occultation light curves produce a data product in 
the form of number of photon counts (flux) vs. time (see 
Elliot & Olkin 1996, for a review). To turn these data into 
a form that is more physically intuitive, the data are fit or 
inverted to obtain temperature vs. altitude (or radius from 
Pluto’s center). This inversion technique can induce large er¬ 
rors between the inverted temperature profile and the true 
temperature profile (as tested with synthetic light curves). 
Model fits assume a certain structure. Thus, when compar¬ 
ing light curve models, it is best to take a temperature profile 
and carry out the forward problem—calculate a model light 
curve and compare data and model in the flux vs. time (or 
flux vs. altitude) regime. 

Chamberlain & Elliot (1997) and Zalucha et al. (2011a) 
describe in detail how to calculate a model light curve. 
Here the procedure is briefly summarized. Given a temper¬ 
ature profile as a function of pressure (in this Case from the 
PGCM), we may use the ideal gas law and the law of hy¬ 
drostatic balance to obtain number density as a function of 
distance from Pluto’s center. The number density is propor¬ 
tional to the refractivity, and using the laws of geometric 
optics, we may obtain the bending angle and its derivative 
with radius. Once these quantities are known, along with 


the distance from the observer to the occulting body, a light 
curve may be calculated. 

Figure 8 shows an example of what ground-based ob¬ 
servers would observe in the year 2015 given different at¬ 
mospheric configurations. Extinction by aerosols or other 
particles has not been taken into account, because including 
them would be inconsistent with the PGCM not accounting 
for their radiative effects; thus, an upturn (i.e., the central 
flash due to diffraction around Pluto’s disk) is present near 
the middle of the light curve. The light curve has been cast 
as intensity vs. altitude in Pluto’s atmosphere. To convert 
a particular observing geometry and the quantities that are 
measured, i.e. flux vs. time, one uses the equations: 

s{r) = r + Dd{r) (14) 

i - (15) 

where s is the observer’s position in the shadow (observer’s) 
plane, r = z + rs (where rs is the body’s radius) is the 
distance from the body’s center in the body (Pluto) plane, 
D is the distance between the observer and the body, 6 
is the bending angle, tmid is the midtime of the occulta¬ 
tion, Smid is the distance of closest approach (i.e., impact 
parameter) in the shadow plane, and V is the relative ve¬ 
locity between the observer and the body. In Fig. 8, there 
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Figure 7. Instantaneous vertical wind (10“®/jbar s“^) for Case 3 on 2015 July 14. Positive values are downward. On this date at the 
time of the New Horizons flyby, the subsolar longitude was 32° and the subsolar latitude was 38° (in the pole convention used in this 
work, see text). Top left panel: 2.2 km altitude, top right panel: 32 km altitude, bottom left panel: 98 km altitude, bottom right panel: 
217 km altitude. On this date, the subsolar point is located at 32° longitude. Except at the lowest level, air appears to be flowing over 
the poles along the terminators. The magnitude of the maximum flow is about 10m s"'^. 




is a clear difference between the light curves derived from 
PGCM simulations with different CH 4 concentrations and 
a slight difference between light curves derived from sim¬ 
ulations with different initial surface pressures. The shape 
of a light curve depends on the temperature gradient with 
height in the atmosphere (Elliot & Young 1992). The be¬ 
haviour of the light curves is consistent with the PGCM 
temperature results (Fig. 2 ), where there is a larger differ¬ 
ence between the cases with different CH 4 concentrations 
than there is between different initial surface pressures. This 
behaviour of the light curves is opposite that of Zalucha et al. 
( 2011 a), who found that when using the steady state version 
of the Strobel et al. (1996) model, the greatest difference ap¬ 
peared in light curves of different surface pressure and the 
difference between CH 4 concentrations was secondary. 

There is some difference in how the Strobel et al. (1996) 
was applied in Zalucha et al. (2011a) and this work. First, 
the original Strobel et al. (1996) version assumed disk av¬ 
eraged radiation, while in the process of implementing it 
into the PGCM, dependence on insolation in the column 
as a function of latitude, longitude, and season was added. 
Second, when the temperature is calculated in the PGCM, 
it is subject to not only the heating and cooling terms of 
Strobel et al. (1996), but also other temperature modifying 
terms: advection, adiabatic heating and cooling, and numer¬ 


ical filters. These in turn interplay with the entire system of 
equations that comprise the PGCM. It is certain that the 
Strobel et al. (1996) model has been correctly implemented 
into the PGCM because if the calls to the dynamic sub¬ 
routines and other subroutines are disabled, the results of 
the stand-alone Strobel et al. (1996) model are recovered. 
In summary, a ground-based observer would be able to dis¬ 
tinguish between CH 4 concentrations of 0.2% and 1% using 
the PGCM, but probably not between pressures of 8 and 
24 pbar. 


6 PREDICTIONS FOR NEW HORIZONS 

When a radio occultation experiment is performed, radio 
waves of a known frequency are transmitted through an at¬ 
mosphere and are measured by a receiver. The rays bend 
as they travel through the atmosphere due to diffraction of 
light in a medium. Unlike shorter wavelengths of light, the 
differential refraction (i.e. flux difference) of radio waves due 
to the density gradient in an atmosphere is too small to be 
measured. What is measured is the difference in Doppler- 
shifted frequencies due to light travelling through the at¬ 
mosphere compared with empty space. The bending angle 
can be calculated given precise knowledge of the geometry 


MNRAS 000, 1-18 (2015) 
















14 


A. M. Zalucha 



Altitude (km) 


Figure 8. Model stellar occultation light curves. Red: initial 
surface pressure 8 fihar and CH 4 concentration 1%, green: ini¬ 
tial surface pressure 24 /rbar and CH 4 concentration 1%, dark 
blue: initial surface pressure 8 /rbar and CH 4 concentration 0.2%, 
magenta: initial surface pressure 24 /rbar and CH 4 concentration 
0.2%. The red and green curves are almost indistinguishable as 
are the blue and magenta curves. A ground-based observer could 
distinguish between CH 4 concentrations of 0.2% and 1% using the 
PGCM, but probably not between pressures of 8 and 24 /rbar. 


of the transmitter and receiver (Hinson et al. 1999) in the 
limit of geometric optics. The bending angle as a function 
of impact parameter relative to the occulting body’s cen¬ 
ter may then be transformed into the refractive index as a 
function of radius from the occulting body’s center. The in¬ 
dex of refraction, rather the refractivity, is proportional to 
the number density, and using the ideal gas law and the as¬ 
sumption of hydrostatic balance, temperature as a function 
of height may be obtained. The solar and stellar occultation 
experiment performed by Alice use the same equations as in 
Section 5. 

Table 3 shows the longitudes and latitudes that each of 
the New Horizons occultation experiments probed. Figs. 9- 
11 show the temperature profiles predicted to be derived 
from by Alice (sun), Alice (background star), and REX ob¬ 
servations. The Alice and Rex profile altitudes were calcu¬ 
lated with the temperature and pressure of that column. 
There is little difference between the immersion and emer¬ 
sion profiles and the different experiments, owing to the weak 
horizontal temperature gradients described in Section 4.1. 
The effect of CH 4 concentration may be seen (c.f. solid 
and dotted, dashed and dot-dashed lined) and initial surface 
pressure (c.f. solid and dashed, and dotted and dot-dashed). 
Near 25 km altitude, the temperature profiles differ by as 
much as 20 K, which will be distinguishable by New Hori¬ 
zons. Above 50 km altitude, the maximum difference is as 
much as 10 K, which is still distinguishable, but in some 
places the difference is zero. However, as we are comparing 
the whole temperature prohle and not individual points, the 
PGCM should be able to tell the difference between differ¬ 
ent initial surface pressures and CH 4 concentrations. Further 
simulations may be needed to pinpoint the combination of 
initial surface pressure and CH 4 concentration, which can 
be carried out if necessary. 


Table 3. Longitude and latitude coordinates for New Horizons 
occultations 


Coordinate 

Alice 

(Sun) 

Alice 

(Background star) 

REX 

Immersion longitude 

150° 

170° 

164° 

Immersion latitude 

0 

0 

0 

0 

15° 

Emersion longitude 

350° 

340° 

343° 

Emersion latitude 

- 20 ° 

- 10 ° 

-14° 



Figure 9. Temperature profiles predicted to be derived from 
Alice observations during the solar occultation. Red: immer¬ 
sion, blue: emersion. Solid line: initial surface pressure 8 fihar 
and CH 4 concentration 1%, dashed line: initial surface pressure 
24 /jbar and CH 4 concentration 1%, dotted line: initial surface 
pressure 8 ^bar and CH 4 concentration 0.2%, dot-dashed line: 
initial surface pressure 24 /rbar and CH 4 concentration 0.2%. The 
immersion and emersion profiles do not differ greatly, which is to 
be expected since the horizontal temperature gradients are weak. 
The temperature profiles for different initial surface pressures and 
CH 4 concentrations differ by up to 20 K in some places, which 
will be detectable by New Horizons. 



Figure 10. Same as Fig. 9, but for the Alice stellar occultation. 
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Altitude (km) 

Figure 11. Same as Fig. 9, but for REX the radio occultation. 

7 DISCUSSION 

Toigo et al. (2015) performed simultaneous, independent 
simulations of Pluto’s atmosphere from this work. They used 
a different approach, by which they ran the surface model 
to equilibrium before coupling it to a GCM. The dynam¬ 
ical core is that of planetWRF (Richardson et al. 2007). 
They used the Zhn et al. (2014) scheme (the snccessor to 
Strobel et al. (1996), mainly adjusted to predict conditions 
at very high altitudes) for atmospheric heating and cooling. 
There are some fundamental differences between the dynam¬ 
ical cores of planetWRF and the MIT GCM. First, plan¬ 
et WRF uses a finite difference scheme while the MIT GCM 
dynamical core uses a hnite volume scheme. Numerically, 
the former approximates derivatives as discrete differences 
between grid points, whereas the latter casts the primitive 
equations of fluid dynamics in terms of fluxes, and calcu¬ 
lates fluxes across grid box boundaries. A second difference 
between planetWRF and the MIT GCM dynamical cores is 
that planetWRF uses a latitude longitude grid in spherical 
coordinates, while the MIT GCM uses a cubed-sphere grid: 
a projection of a cube onto a sphere. The drawback of a 
spherical projection is that near the poles, the grid points 
become very close to each other, eventually turning into a 
numerical singularity at the pole. In order to not violate the 
CFL criterion, a prohibitively small time step is necessary. 
To get around this issue, users of spherical grids use a nu¬ 
merical filter to artificially dampen small scale oscillations. 

The results of this work and Toigo et al. (2015) agree 
in that the latitudinal temperature gradient is weak. How¬ 
ever due to the fact that Toigo et al. (2015) uses much dif¬ 
ferent boundary conditions, which among themselves are 
highly different (following that of Young 2013), they ob¬ 
tain different wind structures. Because at present a cou¬ 
pled GCM and subsurface/surface model has been 
spun-up simultaneously, it cannot be said which initial¬ 
ization method for the lower boundary condition, i.e that of 
this work or those of Toigo et al. (2015), is the proper one 
(though in light of New Horizons Hndings, both will need 
to be modified). It is not known if equilibrating the surface 
frost with a minimally complex atmosphere is adequate or 
if the subsurface/surface and atmosphere are tightly cou¬ 


pled. Nor is it known that if they are tightly coupled, what 
the implications are for the surface frost distribution and 
atmospheric state variables. At this point in time we may 
only describe the effects of the atmosphere given a certain 
boundary condition and compare with observables such as 
temperature and surface frost distribution. This issue will 
only be resolved when a fully coupled PGGM is integrated 
until it has reached equilibrium (determined by repeatabil¬ 
ity of seasonal surface frost and quasi-equilibrium of globally 
averaged atmospheric kinetic energy). 


8 CONCLUSIONS 

This paper has presented results from a PGGM based on the 
MIT GCM. It has continued the work of ZG12 and ZM13, 
who also used the MIT GCM as a basis. This latest version 
uses the radiative-conductive scheme of Strobel et al. (1996) 
for atmospheric heating and cooling, which adds cooling by 
CO and heating by the 2.3 fj,m CH 4 line to the Yelle & 
Lunine (1989) model used in previous versions; mass ex¬ 
change between the surface and atmosphere; a sur¬ 
face and subsurface model; a drag force (to prevent re¬ 
flections off the artificial model top and represent buoyancy 
waves) at the model top that is several orders of magni¬ 
tude weaker than previous versions; and a spin up procedure 
that has Newtonian relaxation step at the beginning to more 
gradually initialize the model, which allows for a larger time 
step (and thus a more efficient use of super-computing time). 
This model of Pluto is one of two very similar models de¬ 
veloped at the same time by independent groups (the other 
being Toigo et al. 2015). 

Importantly, the PGGM can predict wind, a quan¬ 
tity that is difficult to measure remotely even for a satellite 
in orbit or spacecraft flybys. Special circumstances must be 
present, for example relating ocean wave heights to surface 
wind (Earth and Titan), observing tracers such as clouds 
or plumes (Venus, Earth, Triton, and Jupiter), or observing 
wind surface patterns in dunes or ice (Mars and Triton). 
Often these phenomena are transient or only provide wind 
direction and not speed. The fact that the PGGM can pre¬ 
dict winds over the entire globe at any time is a key power 
shared with GCMs in general. 

The model was integrated for 30 Earth-years (0.12- 
Pluto years)from the years 1985 to 2015 to cover the stel¬ 
lar occultation observational record. Due to the computa¬ 
tionally prohibitively expensive spin-up time to equilibrate 
the atmosphere and the surface, it is unknown if the model 
achieved equilibrium. But, based on other modelling and ob¬ 
servational studies it is not likely the case. Thus, four simula¬ 
tions were performed that varied initial global mean surface 
pressure and methane concentration. As in ZM13, there is 
effectively no longitudinal (diurnal) atmospheric temper¬ 
ature variation, and the latitudinal atmospheric temper¬ 
ature gradient is small. The first property follows from the 
estimation that the radiative-conductive time-scale is much 
longer than a Pluto day. However it has yet to be explained 
why a large latitudinal gradient is not predicted on Pluto. 
Given the non-uniform ice distribution and composition ob¬ 
served by New Horizons (Stern 2015), this property is likely 
to change. 

The longitudinal wind has a completely different struc- 
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ture than ZM13 (and ZM12, but this is to be expected since 
that model was 2-D), but the maximum magnitude of order 
10 m s“^ is the same. In the current study, the longitudinal 
winds are of order 1 m s“^ below 10 km. There is a wave 
number -1 structure in this altitude zone in the longitudinal 
direction. Above 10 km altitude upwards to 75 km altitude 
(except in one case where this region extends to 130 km al¬ 
titude), the wind is westward everywhere. Near the south 
pole there is a strong westward jet, representing the global 
maximum wind speed. Above this region up to 175 km alti¬ 
tude, there is a region of weak (~1 m s“^) eastward winds 
between ±30° latitude. Poleward of this latitude the winds 
are westward. From 175 km altitude to the base of the drag 
layer (350 km altitude) there is a wave number-1 structure 
between ±60° latitude with a local maximum of eastward 
winds at 270° longitude and a local maximum of westward 
winds at 90° longitude. Poleward of these latitudes, the lo¬ 
cation of the eastward and westward maxima are shifted by 
180°. 

In the longitudinal average, the latitudinal and verti¬ 
cal winds averaged to approximately zero. However when 
viewed on horizontal planes at different altitudes, it is ap¬ 
parent that air is flowing over the poles (albeit along differ¬ 
ent meridians depending on height), with broad areas of up- 
welling east of the subsolar point and down-welling west of 
the subsolar point. Once again we find Pluto’s atmosphere is 
unlike any other in the solar system. There is a caveat that 
the subsurface/surface and atmosphere may not have had 
time to equilibrate, due to the computationally prohibitive 
wall clock times to run a simulation for even only one Pluto 
year. Thus it is possible the flow will be different, if and 
when it is possible to carry out such long-term simulations. 

Some notable observations were made when simulations 
were repeated with various processes left out. First, when 
the model was initialized at Northern Hemisphere winter 
solstice, the longitudinal winds became eastward almost ev¬ 
erywhere with a maximum above the equator. When the 
contribution from the 2.3 pm CH 4 band was removed, the 
winds became stagnant at the equator. When the upper at¬ 
mosphere drag force was increased and extended to lower 
levels, the top of the longitudinal flow lowered, but the struc¬ 
ture stayed the same, only compressed. Little difference was 
seen when the effects of CO were removed from the radiative- 
conductive scheme. When the surface-atmosphere mass 
exchange was turned off, the longitudinal winds just be¬ 
low the damping layer strengthened at the equator. 

Predictions for the temperature that will be derived 
from New Horizons Alice and Rex observations were pro¬ 
vided. The location (i.e. latitude and longitude) that each 
instrument sampled did not produce significant difference 
in the temperature profiles, since temperature gradients in 
the PGCM results were weak. This includes immersion and 
emersion profiles. In the stratosphere, the temperature pro¬ 
files differ by as much as 20 K, which will be distinguishable 
by New Horizons. In the mesosphere (the more isothermal 
layer above 100 km) the maximum difference is as much as 
10 K, which is still distinguishable, but in some places the 
difference is zero. However, as we are comparing the whole 
temperature profile and not individual points, the PGCM 
should be able to tell the difference between different initial 
surface pressures and CH 4 concentrations. 

In the predictions for ground-based stellar occultations. 


the difference between light curves of different CH 4 con¬ 
centrations would be detectable, but the difference between 
light curves of different initial surface pressures would not. 
While this result is opposite that of Zalucha et al. (2011a), 
we can rule out coding error by de-activating everything but 
the call to the Strobel et al. (1996) module. The difference 
must lie in the fact that the Strobel et al. (1996) module of 
the PGCM has longitudinally, latitudinally, and seasonally 
varying insolation, and that the model is no longer in the 
steady state [dT/dt — 0), but interplay with the rest of the 
model (Navier stokes equations, mass and energy equations, 
ideal gas law, numerical filters). 

Future simulations will test the subsurface and 
surface models and the volatile cycle before conduct¬ 
ing runs spanning multiple Pluto years. It is specu¬ 
lated that this is enough time to spin up the volatile cycle, 
i.e. the location and depth of surface ice, rather than assume 
effectively infinite ice cover. Mars GCMs require about two 
Martian years to equilibrate the ice cycle. From Hansen & 
Paige (1996), it was found that the frost cycle was spun up 
after four Pluto years. This model had a one-layer atmo¬ 
sphere, so it is unknown at this time what the effect of a 
complete surface and atmosphere model will require for spin 
up. A fully spun up volatile cycle might affect the results for 
winds and temperature, particularly the latitudinal wind, 
as air flows from the sublimating summer pole (an area of 
high pressure) to the condensing winter pole (an area of low 
pressure), as is seen in in the lowest atmospheric layer of 
this PGCM and Mars GCMs. 

An interesting insight into volatile transport in the 
PGCM would be to use the “tracer” package. This feature 
allows the user to initialize inert particles in the atmosphere, 
and the scheme keeps track of where they are transported 
(though making these particles radiatively active would not 
be computationally difficult). This treatment would allow us 
to also keep track of sublimated N 2 (including if it froze and 
resublimiated) if say it was transported up to another layer. 
Future work is planned to carry out these simulations. 


ACKNOWLEDGEMENTS 

This work used the Extreme Science and Engineering Dis¬ 
covery Environment (XSEDE), which is supported by Na¬ 
tional Science Foundation Grant No. OGI-1053575. The au¬ 
thor is supported by NASA grant NNX136AH77G. The au¬ 
thor thanks Alan Stern, Dave Hinson, Leslie Young, and 
Henry Throop for providing coordinates of the New Hori¬ 
zons temperature retrievals; Leslie Young for discussions re¬ 
garding the subsurface module; Tim Michaels for discussion 
regarding the computational aspects of the model; Jake Si¬ 
mon for providing comments on the manuscript; and an 
anonymous referee for providing additional comments on the 
manuscript. 

REFERENCES 

Adcroft A., Campin J.-M., 2004, Ocean Modelling, 7, 269 
Adcroft A., Campin J.-M., Hill C., Marshall J., 2004, Mon. Wea. 
Rev., 132, 2845 

Buie M. W., Grundy W. M., Young E. F., Young L. A., Stern 
S. A., 2010, Astron. J., 139, 1117 


MNRAS 000, 1-18 (2015) 


GCM predictions for New Horizons temperatures 17 


Chamberlain D. M., Elliot J. L., 1997, Publ. Astron. Soc. Pacific, 
109, 1170 

Cruikshank D. P., Roush T. L., Moore J. M., Sykes M., Owen 
T. C., Bartholomew M. J., Brown R. H., Tryka K. A., 1997, in 
Stern S. A., Tholen D. J., eds, , Pluto and Charon. University 
of Arizona Press, p. 221 

Elliot J. L., Olkin C. B., 1996, Annu. Rev. of Earth and Planet. 
Sci., 24, 89 

Elliot J. L., Young L. A., 1992, Astron. J., 103, 991 
Elliot J. L., Dunham E. W., Bosh A. S., Slivan S. M., Young 
L. A., Wasserman L. H., Millis R. L., 1989, Icarus, 77, 148 
Elliot J. L., Person M. J., Qu S., 2003a, Astron. J., 126, 1041 
Elliot J. L., et ah, 2003b, Nature, 424, 165 
Elliot J. L., et ah, 2007, Astron. J., 134, 1 
Greaves J. S., Helling C., Friberg P., 2011, MNRAS, 414, L36 
Hansen C. J., Paige D. A., 1996, Icarus, 120, 247 
Held L M., Suarez M. J., 1994, Bull. Amer. Meteor. Soc., 75, 1825 
Hinson D. P., Simpson R. A., Twicken J. D., Tyler G. L., Flasar 
F. M., 1999, J. Geophys. Res., 104, 26997 
Holler B. J., Young L. A., Grundy W. M., Olkin C. B., Cook 
J. C., 2014, Icarus, 243, 104 

Hubbard W. B., Yelle R. V., Lunine J. I., 1990, Icarus, 84, 1 
Hubbard W. B., McCarthy D. W., Kulesa C. A., Benecchi S. D., 
Person M. J., Elliot J. L., Gulbis A. A. S., 2009, Icarus, 204, 
284 

Lellouch E., Sicardy B., de Bergh C., Kaufi H.-U., Kassi S., Cam- 
pargue A., 2009, Astron. Astrophys., 495, L17 
Lellouch E., de Bergh C., Sicardy B., Kaufi H. U., Smette A., 
2011, Astron. Astrophys., 530, L4 
Lewis S. R., Read P. L., Gollins M., 1996, Planet. Space Sci., 44, 
1395 

Marshall J., Adcroft A., Hill C., Perelman L., Heisey C., 1997, J. 
Geophys. Res., 102, 5753 

Michaels T., Young L. A., 2011, in EPSC-DPS Joint Meeting 
2011. p. 1613 

Miller C., Chanover N., Murphy J. R., Zalucha A. M., 2011, AGU 
Fall Meeting Abstracts, p. Al 

Mueller-Wodarg 1. C., Yelle R. V., Mendillo M., Aylward A. D., 
2001, AGU Fall Meeting Abstracts, p. A7 
Nayvelt L., Gierash P. J., Cook K. H., 1997, J. Atmos. Sci., 54, 
986 

Olkin C. B., et ah, 2015, Icarus, 246, 220 
Owen T. C., et ah, 1993, Science, 261, 745 
Pasachoff J. M., et ah, 2005, Astron. J., 129, 1718 
Person M. J., et ah, 2008, Astron. J., 136, 1510 
Reinsch K., Burwitz V., Festou M. C., 1994, Icarus, 108, 209 
Richardson M. L, Toigo A. D., Newman C. E., 2007, J. Geophys. 
Res., 112, 9001 

Showman A. P., Gooper C. S., Fortney J. J., Marley M. S., 2008, 
Astrophys. J., 682, 559 
Sicardy B., et ah, 2003, Nature, 424, 168 

Spencer J. R., Stansberry J. A., Trafton L. M., Young E. F., 
Binzel R. P., Croft S. K., 1997, in Stern S. A., Tholen D. J., 
eds, , Pluto and Charon. University of Arizona Press, pp 435- 
474 

Stansberry J. A., Lunine J. L, Hubbard W. B., Yelle R. V., 
Hunten D. M., 1994, Icarus, 111, 503 
Stern S. A. e. a., 2015, Science, 350, aadl815 
Stern S. A., et ah, 2008, Space Science Reviews, 140, 155 
Strobel D. F., Zhu X., Summers M. E., Stevens M. H., 1996, 
Icarus, 120, 266 

Tholen D. J., Buie M. W., 1990, in Bulletin of the American 
Astronomical Society, p. 7^25.14 

Toigo A. D., French R. G., Gierasch P. J., Guzewich S. D., Zhu 
X., Richardson M. L, 2015, Icarus, 254, 306 
Tryka K. A., Brown R. H., Anicich V., Cruikshank D. P., Owen 
T. C., 1993, Science, 261, 751 


Tryka K. A., Brown R. H., Chruikshank D. P., Owen T. C., 
Geballe T. R., de Bergh C., 1994, Icarus, 112, 513 
Tyler G. L., Linscott 1. R., Bird M. K., Hinson D. P., Strobel 
D. F., Patzold M., Summers M. E., Sivaramakrishnan K., 
2008, Space Science Reviews, 140, 217 
Vangvichith M., Forget F., 2011, in EPSC-DPS Joint Meeting 
2011. p. 1164 

Yelle R. V., Lunine J. L, 1989, Nature, 339, 288 
Young L. A., 2012, Icarus, 221, 80 
Young L. A., 2013, Astro. Phys. J. Letters, 766, L22 
Young E. F., et ah, 2008, Astron. J., 136, 1757 
Zalucha A. M., Gulbis A. A. S., 2012, J. Geophys. Res., 117, 5002 
Zalucha A. M., Michaels T. L, 2013, Icarus, 223, 819 
Zalucha A. M., Gulbis A. A. S., Zhu X., Strobel D. F., Elliot J. L., 
2011a, Icarus, 211, 804 

Zalucha A. M., Zhu X., Gulbis A. A. S., Strobel D. F., Elliot J. L., 
2011b, Icarus, 214, 685 
Zangari A., 2015, Icarus, 246, 93 

Zhu X., Strobel D. F., Erwin J. T., 2014, Icarus, 228, 301 


APPENDIX A: SURFACE AND SUBSURFACE 
ENERGY MODEL 


The subsurface, assumed to be H20(Cruikshank 
et al. (1997) and concluded by Stern (2015), has a 
temperature specified following Young (2012); 

PssCss — Kss 

where psa is the density of the subsurface, Caa is the 
specific heat of the subsurface, Taa is the tempera¬ 
ture of the subsurface, jz is now extended below the 
surface (to negative values), and Kaa is the thermal 
conductivity of the subsurface. A subsurface with 
25 layers is assumed. The thermal inertia T and the 
diurnal skin depth Z are given by 


r — y PasCaaf^aa 


and 


Kis 

ryn 


(A2) 

(A3) 


where Q. — 1.13856 x 10“® s“^ is Pluto’s rotation rate. 
The depth of each layer is 5z = 1/4Z, where the first 
layer begins at 1/2 5z. The depth of the subsurface 
is thus 4.74 m. Using these values, the subsurface 
dynamics captures daily oscillations but not yearly 
ones. However, since these simulations are meant to 
be short-term the impact is assumed to be minimal 
on the final result. Future long-term simulations will 
have a deeper subsurface. 

Equation (Al) is a simple diffusion relationship. 
The discrete form of equation (Al) is 


t: 


rn-|-l 


-T" 


PssCs 


r^+i - 2T7 


+ r"_i 


At 


A 22 


(A4) 


where n indicates the time index and j indicates the 
level. Defining /3 = Kss At/p^sCssAa^, the subsurface 
temperature is time marched by 


T"+i = /3T7+1 + (1 - 2/3) Tf + (A5) 


At the bottom level, it is assumed that no heat is 
being conducted from below, in which case the third 
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term on the right hand side disappears and the “2” 
in the second term on the right hand side becomes 
“1.” Also, for the top level (j = 1), is taken to 

be the snrface temperature. 

The snrface temperatnre, in discrete form, is 
given by 

^ rr + -|^ [s(l - A,) - ^ {T: - TD] 

(A6) 

where Aa is the surface albedo and and pa are the 
specific heat and density of the snrface layer, which 
may either be seasonal N 2 frost or the bare subsur¬ 
face. Note that the subsurface does not snblimate 
or freeze in this model, only surface N 2 . The terms 
in equation (A6) from left to right are insolation S, 
emission by the surface with being the emissivity 
of the surface and a is the Stefan-Boltzmann con¬ 
stant, and conduction between the surface and the 
npper subsurface layer (where Ha is the conductivity 
of the surface). The layer thickness Az, is one fourth 
the skin depth Z = Kss/(rss\/f2). Here Fss is the ther¬ 
mal inertia of the subsurface, equal to sqrtpaaCaai^aa- 

During the 30 Earth-year (0.12 Pluto-year) 
integration time, it was found that the surface 
and snbsurface parameters did not greatly affect 
the results. There is reason to be sceptical of 
this result; multi-Pluto year simulations are in 
progress. For now, the snrface and subsnrface pa¬ 
rameters are as follows: subsurface density and 
subsnrface specific heat capacity (H 2 O ice) are 
936 kg m“® and 960 J kg“^ K“^, respectively, and 
for the seasonal N2(/3) ice, the density, and specific 
heat 1000 kg m~^ and 1340 J kg“^ K“^, respec¬ 
tively (Spencer et al. 1997)^. The condnctivity of 
the snbsurface is 4.88 J m“^ s~^. Until such 

time that multi-Pluto year simulations can be per¬ 
formed and the surface frost distribution can be 
“spun-up” alongside a full atmospheric treatment, 
it is assumed, based on current understanding of 
Pluto at the time of these simulations, that thick¬ 
ness of the N 2 snrface ice is 1 m thick. This value 
was chosen such that there is sufficient N 2 surface 
ice such that it will not sublimate over the integra¬ 
tion time but that the properties of the subsurface 
are still important terms in the surface heat equa¬ 
tion. The surface and subsurface temperature were 
initialized at the N 2 frost temperatnre and isother¬ 
mal with depth. The snrface emissivity was assumed 
to be 1, the surface conductivity 0.2 J m“^ s~^, 

and the surface albedo 0.2. The parameter space ex¬ 
plored by this survey was much larger, a set of cases 
with a surface albedo of 0.9 and another set of cases 
with a surface emissivity of 0.2. All other parameters 


being equal, the results puzzlingly were practically 
independent of surface albedo and emissivity. 

In the GCM results, the surface temperature re¬ 
mained at the N 2 freezing temperature, globally for 
all cases. Over the course of every simulation, N 2 ice 
was found to accumulate everywhere on the surface. 
No gridpoints sublimated the initial surface N 2 ice 
at any time during the simulations. We expect a con¬ 
densation flow (i.e. transport of volatiles by wind), 
but such a flow was not identified. It is likely that the 
surface needs far more than 30 Earth years to equi¬ 
librate, as shown by other modelling studies (e.g. 
Hansen &: Paige 1996; Young 2013). The thickness 
of ice that accumulated was of order 1 mm. No struc¬ 
ture was found in the ice accnmulation vs. latitnde 
and longitnde. The subsurface temperature changed 
by less than 0.01 K at any depth globally. Simula¬ 
tions in progress with integration times of several 
Pluto years will be the subject of future work. More¬ 
over, the subsurface and surface model and the be¬ 
haviour of the volatile scheme will be tested with 
a simplified atmosphere to confirm agreement with 
surface energy balance models (e.g. Hansen & Paige 
1996; Yonng 2013). 


APPENDIX B: ADDITIONAL FIGURES AND 
DISCUSSION FOR CASES 1, 2, AND 4 

Simulations with the the global mean surface pressures and 
CH 4 mixing ratios for Cases 1, 2, and 4 are worth exam¬ 
ining in the came close detail as Case 3 in the main part 
of the results and discussion, because they potentially rep¬ 
resent states that Pluto’s atmosphere may have existed in 
even as recently as the current observational epoch (from 
1988 onwards). Figures B1-B3 show results from the sim¬ 
ulations in which one featnre of the model (CO cooling, 
2.3 pm CH 4 band, upper atmosphere drag coefficient, and 
volatile cycle) was changed, as in Fig. 4. Figures B4-B6 show 
latitude-longitude plots of the longitudinal wind at differ¬ 
ent atmospheric heights, as in Fig. 5. Figures B7-B9 and 
Figs. B10-B12 show latitude-longitude plots of the latitudi¬ 
nal and vertical wind, respectively, at different atmospheric 
heights, as in Figs. 6 and 7, respectively. 

This paper has been typeset from a TpX/UTj;X file prepared by 
the author. 


® A value of 960 J kg“^ was erroneously used for the specific 
heat of H 2 O ice. The value reported by Spencer et al. (1997) is 
350 J kg“^ for ice at 40 K. Given the insensitivity of other 
surface parameters to subsurface temperature discussed later, this 
value is not thought to significantly impact the results; however 
the Spencer et al. (1997) will be used in future work. 
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Figure Bl. Same as Fig. 4 but for Case 1. The panel with the adjacent colorbar is the same as Fig. 3, top left panel for comparison. 
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Figure B2. Same as Fig. 4 but for Case 2. The panel with the adjacent colorbar is the same as Fig. 3, top right panel for comparison. 
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Figure B3. Same as Fig. 4 but for Case 4. The panel with the adjacent colorbar is the same as Fig. 3, bottom right panel for comparison. 
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Figure B4. Same as Fig. 5, but for Case 1. The altitudes of the panels from left to right and top to bottom are 2.4, 37, 109, and 
248 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure B5. Same as Fig. 5, but for Case 2. The altitudes of the panels from left to right and top to bottom are 2.5, 38, 112, and 
249 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure B6. Same as Fig. 5, but for Case 4. The altitudes of the panels from left to right and top to bottom are 2.4, 35, 105, and 
231 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure B7. Same as Fig. 6, but for Case 1. The altitudes of the panels from left to right and top to bottom are 2.4, 37, 109, and 
248 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure B8. Same as Fig. 6, but for Case 2. The altitudes of the panels from left to right and top to bottom are 2.5, 38, 112, and 
249 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure B9. Same as Fig. 6, but for Case 4. The altitudes of the panels from left to right and top to bottom are 2.4, 35, 105, and 
231 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure BIO. Same as Fig. 7, but for Case 1. The altitudes of the panels from left to right and top to bottom are 2.4, 37, 109, and 
248 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure Bll. Same as Fig. 7, but for Case 2. The altitudes of the panels from left to right and top to bottom are 2.5, 38, 112, and 
249 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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Figure B12. Same as Fig. BIO, but for Case 4. The altitudes of the panels from left to right and top to bottom are 2.4, 35, 105, and 
231 km. (The difference from the other simulations results from the conversion from pressure to altitude). 
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